The Ontogeny and Evolution of Sex-Biased Gene Expression in 
Drosophila melanogaster 

Jennifer C Perry/'^'^'^ Peter W. Harrison,^ and Judith E. Mank^ 

^]esus College, University of Oxford, Oxford, United Kingdom 

^Edward Grey Institute, Department of Zoology, University of Oxford, Oxford, United Kingdom 
^Department of Genetics, Evolution and Environment, University College London, London, United Kingdom 
*Corresponding author: E-mail: jennifer.perry@zoo.ox.acuk. 
Associate editor: Jianzhi Zhang 



Abstract 

Sexually dimorphic phenotypes are thought to largely result from sex differences in gene expression, and genes with sex- 
biased expression have been well characterized in adults of many species. Although most sexual dimorphisms manifest in 
adults, many result from sex-specific developmental trajectories, implying that juveniles may exhibit significant levels of 
sex'biased expression. However, it is unclear how much sex-biased expression occurs before reproductive maturity and 
whether preadult sex-biased genes should exhibit the same evolutionary dynamics observed for adult sex-biased genes. In 
order to understand the continuity, or lack thereof, and evolutionary dynamics of sex-biased expression throughout the 
life cycle, we examined sex-biased genes in pre-gonad tissue of two preadult stages and compared them with the adult 
gonad of Drosophila melanogaster. We found that the majority of the genome is sex-biased at some point in the life cycle, 
with some genes exhibiting conserved sex-biased expression and others displaying stage-specific sex bias. Our results also 
reveal a far more complex pattern of evolution for sex-biased genes throughout development. The most rapid evolu- 
tionary divergence occurred in genes expressed only in larvae within each sex, compared with continuously expressed 
genes. In females — but not males — this pattern appeared to be due to relaxed purifying selection in larva-limited genes. 
Furthermore, genes that retained male bias throughout life evolved more rapidly than stage-specific male-biased genes, 
due to stronger purifying selection in stage-specific genes. However, female-biased genes that were specific to larvae 
evolved most rapidly, a pattern that could not be definitively attributed to differences in adaptive evolution or purifying 
selection, suggesting that pleiotropic constraints on protein-coding sequences can arise when genes are broadly expressed 
across developmental stages. These results indicate that the signature of sex-specific selection can be detected well before 
reproductive maturity and is strongest during development. 
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Introduction 

Many sexual species exhibit considerable sexual dimorphism. 
These sex differences are often complex and therefore involve 
many genes. However, males and females share the majority 
of their genome, suggesting that much of the basis for this 
dimorphism must result from differences in expression of the 
same genes. Sex differences in gene expression have been well 
characterized in adult animals, particularly in classic model 
organisms, and several patterns have emerged that appear 
common to a wide range of taxa. In adults, more genes show a 
stronger bias toward male expression than toward female 
expression (reviewed by Ellegren and Parsch 2007; Zhang 
et al. 2007; Small et al. 2009; Graveley et al. 2011; Zhao et al. 
2011; Martins et al. 2013). Furthermore, male-biased genes 
typically show more rapid gene turnover compared with 
female-biased and unbiased genes, and elevated sequence 
divergence, due at least in part to more rapid adaptive evo- 
lution, particularly for genes expressed in the germline. In 
contrast, female-biased genes often evolve at similar or 
slower rates compared with unbiased genes (Ellegren and 
Parsch 2007; Zhang et al. 2007; Parsch and Ellegren 2013). 



Together, these observations have led to the hypothesis 
that much of the evolution of sex-biased genes is driven by 
selection on males via sexual selection and conflict (Swanson 
and Vacquier 2002; Ranz et al. 2003; Connallon and Knowles 
2005; Small et al. 2009; Artieri and Singh 2010). 

Currently, studies of sex-biased gene expression focus 
almost exclusively on adults, and the ontogenetic complexity 
of sexual dimorphism and sex-biased genes remains largely 
unexplored. This trend reflects the expectation that sex- 
biased expression and sex-specific selection will be minimal 
in juveniles, which is in many ways reasonable, based on the 
outward apparent sexual monomorphism in juveniles of 
many species. Moreover, empirical studies suggest that the 
sexes are subject to similar selection as juveniles, with sex- 
specific selection and intralocus conflict only beginning in 
adults when the sexes have divergent reproductive interests 
(Chippindale et al. 2001; Rice and Chippindale 2001; Gibson 
et al. 2002). Thus, to the extent that sex-biased gene expres- 
sion represents a history of sexual antagonism, minimal juve- 
nile sex-biased expression should be expected, and those 
genes that show sex-biased expression only in juveniles 
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should not display the same rapid evolution exhibited by 
adult sex-biased genes. 

However, recent studies have challenged the generality 
both of juvenile sexually monomorphic gene expression 
and of uniform selection on juvenile males and females. 
A few studies have reported considerable sex-biased gene 
expression in juvenile stages (e.g., rainbow trout. Hale et al. 
2011; silkworm, Zhao et al. 2011), and in some cases the pat- 
tern of sex-biased expression is quite distinct from the adult 
stage (Mank et al. 2010). Furthermore, sex-specific selection 
can act on juvenile fitness even without gross phenotypic 
differences between the sexes (Zikovitz and Agrawal 2013). 
These results are hard to reconcile with the view of juvenile 
sexual monomorphism. As a result, we currently lack a com- 
prehensive understanding of the patterns and evolutionary 
dynamics of sexually dimorphic gene expression throughout 
development. 

Here, we use high-throughput RNA sequencing to investi- 
gate the ontogenetic basis of sex-biased gene expression in 
Dmsophila melanogaster and ask whether the expression and 
evolutionary dynamics of sex-biased genes in two preadult 
stages (third instar larvae and pre-pupae) reflect the adult 
pattern. In D. melanogaster, male-biased genes show a greater 
magnitude of bias than female-biased genes (Ranz et al. 2003; 
Ellegren and Parsch 2007) and exhibit more rapid protein 
divergence, due at least in part to more rapid adaptive evo- 
lution (Zhang et al. 2004; jagadeeshan and Singh 2005; Zhang 
and Parsch 2005; Gnad and Parsch 2006; Proschel et al. 2006; 
Haerty et al. 2007; Zhang et al. 2007; Baines et al. 2008; Artieri 
et al. 2009; Meisel 2011; Grath and Parsch 2012; Muller et al. 

2012) . These patterns are related to both sex differences in ex- 
pression breadth across tissues, with sex-biased genes specific 
to reproductive tissues evolving more rapidly than broadly 
expressed genes (Grath and Parsch 2012; Parsch and Ellegren 

2013) , and chromosomal linkage, with male-biased X-linked 
genes evolving most rapidly (Baines et al. 2008; Grath and 
Parsch 2012; Kayserili et al. 2012; Muller et al. 2012; Parsch and 
Ellegren 2013; see also Meisel et al. 2012b). 

To examine sex-biased expression, we focus on the most 
sexually dimorphic tissue of pre-gonad imaginal discs (in 
larvae and pre-pupae, stages at which sex differentiation is 
an ongoing process) and adult gonads. By isolating a single 
tissue, our study overcomes the problem of tissue allometry 
that has complicated analyses of sex-biased gene expression 
in studies of whole organisms, which can create an imprecise 
picture of overall sex bias (Chintapalli et al. 2007; Catalan et al. 
2012; Perry and Mank forthcoming) and which very likely 
varies across larval, pre-pupal, and adult body plans. We ask 
three questions. 1) What is the pattern of sex-biased gene 
expression in pre-adults, including linkage, tissue specificity, 
and gene function? 2) What is the pattern of continuity of 
expression throughout development for male- and female- 
biased genes? 3) How does ontogeny interact with sex- 
specific selection to shape the divergence of sex-biased 
genes? Our results indicate that sex-biased gene expression 
in juvenile stages is both considerable and dynamic, and that 
the evolution of sex-biased genes throughout ontogeny is 
considerably more complex than previously realized. 



Results 

At both the third instar larval and the pre-pupal stages we 
sampled, pre-gonad tissue is more differentiated in males 
than females. In males, primary spermatocytes begin to dom- 
inate the pre-gonads of third instar larvae and are the dom- 
inant cell type in pre-pupae, with the first meiotic divisions 
occurring around the onset of pupation (Bodenstein 1950). 
Spermatocytes synthesize the RNA required for later sperma- 
tid development (Lindsley and Tokuyasu 1980). In females, 
pre-gonad tissue contains only oogonia with no oocytes yet 
present; the first meiotic divisions occur as the adult female 
emerges (King 1970; Mahowald and Kambysellis 1980). 
Ovarioles are well differentiated in pre-pupae, but differenti- 
ation has not yet commenced at the larval stage we sampled 
(Bodenstein 1950). 

We found that gene expression measures among our bio- 
logical replicates of each sex and life stage were strongly and 
positively correlated, with r ranging from 0.947 to 0.999 
(P < 0.00005). Adult gene expression was positively correlated 
with a previous adult data set (r/v\ALE = 0.865, n - 8,847, 
P < 0.0001; rpEMALE = 0.846, n = 7,474, P < 0.0001; Can et al. 
2010). 

Sex-Biased Gene Expression in Larvae and Pre-pupae 
In contrast to previous whole-body studies of pre-adult 
Dmsophila (Lebo et al. 2009; Meisel et al. 2012a), we found 
that more than half of all expressed genes showed moderate 
to high levels of sex-biased expression (>2-fold expression 
difference, Padj < 0.05). In both larvae and pre-pupae, more 
genes were female- than male-biased for all but the most 
extremely sex-biased genes (those with > 10-fold expression 
difference, Padj < 0.05; fig. 1 and supplementary table SI, 
Supplementary Material online). However, the magnitude 
of differential expression between the sexes was greatest for 
male-biased genes (fig. 1 and supplementary table SI, 
Supplementary Material online). These patterns were broadly 
similar in larvae and pre-pupae. 

Compared with autosomes, a larger proportion of genes 
on the X chromosome was female-biased and a lower pro- 
portion was male-biased in both larvae and pre-pupae (fig. 1 
and supplementary table SI, Supplementary Material online), 
consistent with both the feminization of the X chromosome 
(Meisel et al. 2012a) and incomplete or absent dosage com- 
pensation in the male germline (Meiklejohn et al. 2011; 
Meiklejohn and Presgraves 2012; but see also Hense et al. 
2007; Vicoso and Charlesworth 2009). To further assess 
whether dosage compensation was absent in our preadult 
samples, we calculated 1) the ratio of male to female expres- 
sion for X-linked genes at each stage (log2(mean male expres- 
sion)/log2(mean female expression)), which should equal one 
with full dosage compensation, and 2) the ratio of X to 
autosomal expression in each sex and at each stage, which 
should equal one with full dosage compensation. We found, 
first, that the malefemale expression ratio was 0.85 in larvae 
and 0.79 in pre-pupae. Second, the male X:autosome ratios 
were 0.90 and 0.88 in larvae and pre-pupae, respectively. 
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Fic. 1. Relative gene expression in males and females, for larval and pre-pupal autosomal and X-linked genes. Percentages within each panel indicate the 
percentage of all genes expressed at that stage that were at least 2-fold male-biased (MB; upper left, top number), male-limited (upper left, bottom 
number; i.e., genes with no expression in females), female-biased (FB; lower right, top number), and female-limited (lower right, bottom number), with 
all categories further defined by significantly different expression between the sexes (Pa^j < 0.05, FDR-corrected), The number of expressed genes is given 
in the lower left. Unbiased (UB) genes were defined as Padj > 0.05 and less than 2-fold difference between the sexes. 



The data therefore support a scenario of incomplete dosage 
compensation at these developmental stages. 

Extreme sex bias (>10'fold difference, Padj < 0.05) arises in 
different ways for male- and female-biased genes (fig. 2). 
Extremely male-biased genes are the product of both increased 
expression in males and decreased expression in females. In 
contrast, extremely female-biased genes are not upregulated in 
females but instead arise from down regulation in males. This 
pattern was more pronounced in pre-pupae than in larvae. 

Continuity of Sex-Biased Gene Expression 
Relative gene expression was highly conserved between larvae 
and pre-pupae (supplementary fig. SI, Supplementary 
Material online). Within a sex, the developmental stages 
showed much greater similarity in gene expression than 
the sexes showed with each other within a developmental 

stage (/^[FEMALE LARVAE-PREPUPAE] = 0.8831; i^[MALE LARVAE- 

PREPUPAE] = 0.9553; r [larval male-female] = 0.3984; r [prepupal 
MALE-FEMALE] = 0.2737; all P < 0.0001). 



To examine sex-bias conservation through to the adult 
stage, we examined the overlap in sex-bias categories 
among stages (fig. 3 and supplementary figs. SI and S2, 
Supplementary Material online). For autosomal loci, male- 
biased genes more often retained their sex bias throughout 
life (58% vs. 46% for female-biased genes), whereas female- 
biased genes more often showed stage-specific sex bias (30% 
vs. 25%; fig. 3). This pattern was reversed for X-linked loci, 
potentially due to the overall reduction in male-biased genes 
on the D. melanogaster X. Adults showed higher stage speci- 
ficity for sex bias compared with larvae and pupae. 
Altogether, 46-65% of sex-biased genes retained their sex 
bias throughout life, whereas 18-30% showed stage-specific 
sex bias. 

For the above analysis, we categorized adult sex bias using 
our adult data (see Materials and Methods). We compared 
the results with those obtained when adult sex bias was cat- 
egorized using a meta-analysis of previous gene expression 
studies on mainly whole body flies (SEBIDA; Gnad and Parsch 
2006; supplementary fig. S2, Supplementary Material online). 
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Fig. 2. Mean relative gene expression (fragments per kilobase of transcript per million mapped reads) at different thresholds of sex-bias, for male- and 
female-biased genes expressed in males (blue) and females (red). For each threshold, sex-biased genes were categorized by a significant difference 
between the sexes (FDR-corrected). Error bars were too small to be meaningfully displayed. 



We focused on the meta-analysis data as the best available 
summary across many studies. Previous individual studies 
that have examined adult gonad alone (as we have; e.g., 
Parisi et al. 2003, 2004) have been based on microarrays 
with less sensitivity to lowly expressed genes (Harrison, 
Wright et al. 2012), therefore resulting in relatively few 
genes called as sex-biased, especially for female-biased 
genes. In fact, by comparing data from five studies of adult 



D. melanogaster (including this one), we found that gene 
expression data are more similar between whole-body and 
gonad-only studies than between microarray and RNA-seq 
platforms (supplementary table 52, Supplementary Material 
online). Although our data were strongly correlated with 
those of the SEBIDA meta-analysis, our data yielded more 
genes categorized as sex-biased in adults, and therefore 
more genes with overlapping sex-biased expression among 
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Stages (compare fig. 3 and supplementary fig. SI, 
Supplementary Material online). This difference might 
result from the fact that the meta-analysis incorporates 
many studies with multiple sources of variation, or from 
our data being exclusively from the gonad, as some of the 
patterns of sex bias from the gonad can be diluted through 
the inclusion of somatic tissue in whole-body studies, or from 
differences in platforms among studies (Harrison, Wright 
et al. 2012). 

Reversals in Sex-Specific Gene Expression 
In order to identify genes that reversed the direction of sex- 
specific expression between developmental stages, we filtered 
our data for genes that were sex-biased (>2-fold differential 
expression between the sexes and Padj < 0.05) or sex-limited 
(expressed in one sex only; P^dj < 0.05) in larvae or pre-pupae, 
and sex-biased or sex-limited in the opposite direction in 
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c-biased gene expression in larvae, pre-pupae, and adults. 

adults. Autosomal genes were more than three times as 
likely to switch from female-biased or female-limited expres- 
sion to male-biased or male-limited expression as vice versa 
(7.5% [182/2,427] vs. 1.8% [37/2,026]; x] = 80.8, P < 0.000001; 
supplementary fig. S3, Supplementary Material online). 
X-linked genes showed a similar pattern, although the differ- 
ence was smaller and not statistically significant (3.3% [25/ 
743] vs. 1.5% [5/333]; x] = 2.9, P = 0.08; supplementary fig. S4, 
Supplementary Material online). For these results, we classi- 
fied adult sex bias using our own data, and we obtained 
similar results using the SEBIDA meta-analysis (supplemen- 
tary text. Supplementary Material online; Gnad and Parsch 
2006). 

Evolutionary Divergence of Sex-Biased Genes 

To test for differences in the rate of evolutionary divergence 

among categories of sex-biased and sex-limited genes, we 
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Table 1. Divergence Estimates for Sex-Biased and Sex-Limited Gene Categories in 


Larvae and Pre-pupae. 






Stage Location 


Category^ 


n Genes (kb/1000)'' 




^^s {PT 


dulds {Pf 


DoS {Pf 


Larva Autosomes 


UB 


659 (41.49) 


0.0062 


0.0662 


0.0930 


-0.14 




MB 


1202 (90.13) 


0.0082 (<0.007) 


0.0744 (<0.007) 


0.1107 {0.004) 


-0.10 (<0.0007) 




ML 


290 (24.83) 


0.0092 (<0.007) 


0.0780 (<0.007) 


0.1179 {0.004) 


-0.12 (0.027) 




FB 


1915 (167.31) 


0.0051 (0.007) 


0.0691 (0.073) 


0.0742 (<0.007) 


-0.16 (0.116) 




FL 


179 (12.23) 


0.0067 (0.410) 


0.0711 (0.023) 


0.0944 (0.892) 


-0.17 (0.113) 


X 


UB 


58 (4.42) 


0.0091 


0.0809 


0.1123 


-0.10 




MB 


111 (11.70) 


0.0149 (0.076) 


0.0975 (0.058) 


0.1524 (0.101) 


-0.06 (0.244) 




ML 


35 (3.85) 


0.0140 {0.044) 


0.1162 {0.004) 


0.1207 (0.811) 


—0.09 (0.858) 




FB 


442 (27.38) 


0.0050 (<0.007) 


0.0746 (0.170) 


0.0676 (0.004) 


—0.14 (0.134) 




FL 


36 (2.51) 


0.0056 (0.029) 


0.0794 (0.845) 


0.0703 (0.083) 


—0.14 (0.458) 


Pre-pupa Autosomes 


UB 


694 (44.31) 


0.0065 


0.0674 


0.0971 


-0.15 




MB 


1029 (75.35) 


0.0080 (<0.007) 


0.0738 (<0.007) 


0.1080 (0.076) 


-0.11 (<0.0007) 




ML 


440 (38.11) 


0.0098 (<0.007) 


0.0796 (<0.007) 


0.1225 (0.007) 


-0.10 (<0.0007) 




FB 


2080 (174.51) 


0.0050 (<0.007) 


0.0682 (0.424) 


0.0733 (<0.007) 


-0.16 (0.214) 




FL 


204 (14.39) 


0.0070 (0.455) 


0.0729 {O.OOG) 


0.0958 (0.881) 


-0.17 (0.074) 


X 


UB 


62 (4.94) 


0.0088 


0.0800 


0.1098 


-0.13 




MB 


95 (9.18) 


0.0120 (0.127) 


0.0961 {0.04G) 


0.1246 (0.490) 


-0.10 (0.109) 




ML 


55 (6.54) 


0.0187 (<0.007) 


0.1135 (<0.007) 


0.1651 (0.071) 


-0.03 (0.005) 




FB 


459 (28.83) 


0.0052 (<0.007) 


0.0755 (0.345) 


0.0690 {0.003) 


-0.14 (0.889) 




FL 


30 (2.06) 


0.0048 {0.038) 


0.0762 (0.601) 


0.0633 (0.081) 


-0.10 (0.314) 



^Unbiased (UB), male-biased (MB), female-biased (FB), male-limited (ML), or female-limited (FL). 
''The number of genes and total alignment length. 

values from a permutation test comparing each set of genes with the unbiased genes for that category. Significant P values (<0.05) are italicized. 
■^Mean values for the direction of selection statistic, calculated as the difference in the proportion of fixed differences that are nonsynonymous and the proportion of 
polymorphisms that are nonsynonymous (Stoletzki and Eyre- Walker 2011). Negative values indicate purifying selection; positive values indicate adaptive evolution. P values 
from a Wilcoxon nonparametric test. 



compared the ratio of nonsynonymous to synonymous nu- 
cleotide substitutions (d^/ds) for sex-biased, sex-limited, and 
unbiased genes in larvae and pre-pupae along the D. melano- 
gaster terminal lineage (table 1). The patterns we observed are 
broadly consistent with previous reports for adult flies; how- 
ever, we identified several interesting deviations. Female- 
biased genes evolved more slowly than unbiased genes in 
both life stages (table 1). In contrast, male-biased and male- 
limited genes evolved more rapidly than unbiased genes, al- 
though this difference was not significant for all categories 
(table 1). 

To assess whether differences in divergence rates were due 
to increased positive selection or relaxed purifying selection, 
we used existing data on standing variation in protein-coding 
sequences (from the Drosophila Genome Reference Panel, 
DGRP; Mackay et al. 2012) to calculate the direction of 
selection (DoS) statistic, a measure of the difference in the 
proportions of nonsynonymous fixed differences and poly- 
morphisms (Stoletzki and Eyre- Walker 2011). Positive DoS 
values indicate positive selection and adaptive evolution, 
whereas zero values indicate neutral evolution and negative 
values indicate purifying selection and the presence of segre- 
gating weakly deleterious mutations. Estimates of DoS were 
negative across all sex-bias categories (table 1). Comparisons 
among categories revealed evidence for relaxed purifying se- 
lection driving the rapid diversification of male-biased and 
male-limited autosomal genes (table 1). In contrast, female- 
biased genes, which evolved more slowly than unbiased genes. 



showed a nonsignificant trend toward stronger purifying 
selection compared with unbiased genes (table 1). 

Evolutionary Divergence of Stage-Specific and 
Continuously Sex-Biased Genes 

Given the extensive overlap in sex bias throughout life (fig. 3 
and supplementary fig. SI, Supplementary Material online), 
the divergence of sex-biased genes in one stage may be driven 
by selection acting at another stage. To examine this possi- 
bility, we tested for differences in evolutionary divergence in 
genes that showed stage-specific or continuous sex-specific 
expression. We focused on autosomal genes because there 
were not enough X-linked stage-specific genes for a meaning- 
ful test. 

We conducted two complementary tests. First, we exam- 
ined divergence in genes with stage-specific expression in 
each sex — that is, expressed in a sex in one developmental 
stage and not expressed in that sex in the other stages — 
compared with genes continuously expressed in that sex 
throughout life. In both sexes, genes with stage-specific 
expression showed more rapid divergence than genes with 
continuous expression in that sex (fig. 4a and supplementary 
table S3, Supplementary Material online). Genes with expres- 
sion limited to the larval stage evolved most rapidly within 
both sexes. We compared DoS estimates among these groups 
of genes to assess the contributions of positive and purifying 
selection to their evolutionary divergence. In females, the 
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Fig. 4. Divergence rates (a, b; the ratio of nonsynonymous to synonymous substitutions, df^/ds) and the direction of selection (c, d) for autosomal genes 
with continuous or stage-limited expression within a sex (a, c), or continuous or stage-limited sex-bias in expression (b, d). Asterisks indicate a significant 
difference (P < 0.05, permutation test) between stage-specific genes and continuously expressed or sex-biased genes for that category. Bars indicate 95% 
bootstrapped confidence intervals; the bars for genes always expressed in {a) are present but very small due to the large sample sizes. The number of 
genes in each category is indicated. 



results are consistent with more rapid evolution due to 
relaxed purifying selection acting on larval-specific genes 
(fig. 4c). In contrast, we found similar DoS values for larval- 
specific and continuously expressed genes in males, and we 
were therefore unable to identify whether the more rapid 
evolution of larval-specific genes was due to relaxed purifying 
selection (fig. 4c). 

Second, we examined genes for which sex bias in expres- 
sion was stage-specific, compared with genes that were con- 
sistently male- or female-biased, or unbiased, throughout life. 
Male- and female-biased genes show distinct patterns of 
evolutionary divergence depending on the extent of sex- 
bias conservation throughout development (fig. 4b and sup- 
plementary table S4, Supplementary Material online). Genes 
that remained male-biased throughout life evolved more 



rapidly than genes that were male-biased in one stage only 
(fig. 4b). This pattern was similar to that of unbiased genes. In 
fact, when this ontogenetic pattern of sex bias was consid- 
ered, male-biased genes no longer evolved more rapidly than 
unbiased genes; the d^/ds ratios of continuously male-biased 
and unbiased genes were statistically indistinguishable (per- 
mutation test, P = 0.605). 

In contrast, female-biased genes with stage-specific sex 
bias evolved more rapidly than genes that maintained their 
female bias throughout development (fig. 4b). The most rap- 
idly evolving female-biased genes — those with larva-specific 
female bias — had a d^/d^ ratio that was statistically indistin- 
guishable from the most rapidly evolving male-biased and 
unbiased genes (with continuous male-biased or unbiased 
expression; permutation tests, P = 0.434, P = 0.808, 
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respectively), although the limited number of genes with 
larva-specific female bias would have reduced the statistical 
power of this test. 

The DoS value for continuously male-biased genes is con- 
sistent with the action of relaxed purifying selection promot- 
ing their more rapid evolution (fig. 4d). However, DoS values 
for consistent and stage-specific genes with either unbiased or 
female-biased expression were statistically indistinguishable 
(fig. 4d), providing no support for differences in positive or 
purifying selection as an explanation for variation in the 
divergence rates of these genes. 

We conducted the above analysis by categorizing adult 
sex-biased genes using our data. When we categorized adult 
sex bias using a meta-analysis of mostly whole body samples 
(SEBIDA, Gnad and Parsch 2006), the results showed similar- 
ities and differences (supplementary fig. S3 and table S5, 
Supplementary Material online). Again, the differences may 
be due to the greater sensitivity of RNA-seq compared with 
microarrays and our use of gonad tissue compared with 
mainly whole-body studies. Using the SEBIDA data, male- 
biased and female-biased genes were distinct from each 
other, with genes that were male-biased in all three life 
stages again evolving more rapidly than genes and with 
DoS values supporting relaxed purifying selection as the driv- 
ing force. Genes that were continuously female-biased again 
evolved relatively slowly. However, here adult-specific female- 
biased genes evolved more rapidly than other female-biased 
genes, with DoS values supporting relaxed purifying selection 
driving this divergence, and unbiased genes showed no 
significant differences in divergence or DoS values between 
continuous and stage-specific expression. 

Evolutionary Divergence of Tissue-Limited and 
Broadly Expressed Genes 

In adult D. melanogaster, female-biased genes are more likely 
than male-biased genes to be expressed in multiple tissues 
(Grath and Parsch 2012). This may explain why female-biased 
genes in adults evolve more slowly, as the constraints of pleio- 
tropic expression in multiple tissues can slow down protein 
evolution (Mank et al. 2008; Meisel 2011). Additionally, adult 
sex-biased genes with tissue-limited expression — particularly 
adult male-biased genes expressed only in reproductive tis- 
sues — tend to evolve more rapidly than those expressed 
broadly (Mank et al. 2008; Mank 2009; Meisel 2011; Grath 
and Parsch 2012). We tested whether these patterns are rep- 
licated in larvae by referencing gene expression data for seven 
larval somatic tissues (available in the FlyAtlas database, 
Chintapalli et al. 2007). We categorized sex-biased and sex- 
limited genes based on whether their expression was limited 
to larval pre-gonad only (i.e., expressed in our data set and 
with signal intensity <100 in each somatic tissue) or larval 
pre-gonad and at least one somatic tissue. Overall, 7,021 genes 
expressed in our larval data set were also reported in the 
FlyAtlas database. Compared with unbiased genes, female- 
biased genes were more likely to be expressed in both larval 
soma and pre-gonad, whereas female-limited, male-biased, 
and male-limited genes were less likely to be expressed 



Table 2. Breadth of Expression of Genes Expressed in Larval Pre- 
gonads, by Sex-Bias Category. 

Location Category^ n Genes Expressed in xl ^ Value 
Both Larval Soma and 
Pre-gonad/Total (%) 



Autosomes UB 558/843 (66%) 

MB 728/1638 (44%) 110.6 <0.000001 

ML 78/437 (18%) 275.1 <0.000001 

FB 2070/2583 (80%) 73.9 <0.000001 

FL 128/297 (43%) 53.3 <0.000001 

X UB 55/116 (47%) 

MB 61/247 (25%) 22.2 0.000015 

ML 7/78 (9%) 35.7 <0.000001 

FB 554/717 (77%) 49.6 <0.000001 

FL 25/65 (38%) 2.8 0.24 



^Genes categorized by their expression in larval pre-gonad as unbiased (UB), male- 
biased (MB), female-biased (FB), male-limited (ML), or female-limited (FL). For each 
category, the frequency of broadly expressed genes is tested against unbiased genes 
for that chromosomal location. 



broadly (table 2). Furthernnore, autosonnal genes were more 
likely to be expressed in multiple tissues than X-linked genes 
(61% vs. 57%; xl = 9.5, P = 0.009). 

To test whether expression breadth influences evolution- 
ary divergence for preadult sex-biased genes, we compared 
genes with expression limited to larval pre-gonads or expres- 
sion in pre-gonads and at least one somatic tissue. In larvae, 
autosomal genes specific to pre-gonad tissue evolved more 
rapidly than genes expressed more broadly (fig. Sa). For most 
sex-bias categories, DoS values did not differ significantly de- 
pending on expression breadth (fig. Sb). A curious exception 
was a higher DoS value for female-biased genes expressed in 
multiple tissues compared with pre-gonad specific genes, 
indicating relaxed purifying selection despite their slower di- 
vergence in protein-coding sequence (fig. 5). However, the 
magnitude of difference in DoS was not large between the 
categories. Like autosomal genes, male- and female-biased 
X-linked genes showed more rapid divergence when they 
were broadly expressed rather than limited to pre-gonad 
tissue (supplementary table S6, Supplementary Material 
online). There were too few X-linked unbiased and sex-limited 
genes with divergence data available to perform a meaningful 
test. 

Gene Ontology 

Gene ontology (GO) analyses can be a useful exploratory tool, 
although it is difficult to draw firm conclusions about func- 
tional differences based on GO analyses alone (e.g., Pavlidis 
et al. 2012). To explore how gene function relates to devel- 
opmental patterns of sex-biased expression, we tested for 
enriched GO terms in larvae and pre-pupae for 1) sex- 
biased, sex-limited, and unbiased genes, 2) genes with contin- 
uous or stage-specific sex bias, and 3) genes with consistent or 
reversed sex bias throughout development. We compared 
each subset of genes to all genes expressed at that stage. 

In both larvae and pre-pupae, female-biased and female- 
limited genes were most enriched for metabolic, cellular, and 
developmental processes, as well as channel activity and 
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Fig. 5. Divergence rates (a; the ratio of nonsynonymous to synonymous substitutions, c/n/c/s) and the direction of selection (b) for unbiased (UB), male- 
biased (MB), male-limited (ML), female-biased (FB), and female-limited (FL) autosomal genes expressed in multiple larval tissues (pre-gonad and at least 
one somatic tissue) or pre-gonad only. Asterisks indicate significant differences (P < 0.05, permutation test). Lines indicate 95% bootstrapped confi- 
dence intervals. The number of genes in each category is indicated. 



compound binding functions, whereas male-biased genes 
were most enriched for male gamete generation and sperma- 
togenesis (supplementary table S6, Supplementary Material 
online). This pattern is consistent with the hypothesis that 
female-biased genes function in more varied ways than male- 
biased genes, which is consistent with female-biased genes 
showing greater breadth of expression across tissues and 
developmental stages. 

Genes that maintained their sex bias throughout develop- 
ment showed distinct enrichment from genes with stage- 
specific sex bias (supplementary table S7, Supplementary 
Material online), particularly for female-biased genes. Genes 
with female bias throughout life were most enriched for cel- 
lular and intracellular parts and cellular processes, whereas 
continuously male-biased genes were most enriched for male 
gamete generation and spermatogenesis. Genes that switched 
from female-biased to male-biased expression throughout 
development were most enriched for extracellular compo- 
nents, whereas genes that switched from male-biased to 



female-biased expression were not significantly enriched for 
any terms (supplementary table S8, Supplementary Material 
online). 

Discussion 

This study expands ongoing efforts to characterize the 
D. melanogaster developmental transcriptome (Arbeitman 
et al. 2002; Daines et al. 2011; Graveley et al. 2011; 
Mikhaylova and Nurminsky 201 1) and to investigate the evo- 
lutionary genetics of potentially sexually antagonistic genes 
(Innocenti and Morrow 2010). Characterizing developmental 
patterns of sex differences in gene expression represents a 
major gap in understanding the evolutionary genetics of 
sex-biased genes (Artieri et al. 2009; Artieri and Singh 2010; 
Carreira et al. 2011; Meisel et al. 2012a; Parsch and Ellegren 
2013). Our study redresses this gap by examining patterns of 
sex-biased gene expression and sex-specific selection in two 
juvenile stages of D. melanogaster. We find that investigating 
sex-biased genes from a developmental perspective yields 
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novel insights. Many genes exhibited conserved sex-biased 
expression throughout development, whereas many others 
displayed stage-specific sex-biased expression. The pattern of 
continuity matters for the evolution of these genes: we de- 
tected more rapid evolution of genes expressed only in larvae 
within each sex compared with continuously expressed genes. 
Furthermore, genes with continuously male-biased expression 
evolved more rapidly than stage-specific male-biased genes 
by relaxed purifying selection, whereas female-biased genes 
showed the reverse pattern but with no evidence for differ- 
ences in purifying or positive selection. Overall, when consid- 
ering gene expression throughout development, differences in 
evolutionary divergence rates between male-biased and 
female-biased genes are minimized. 

Sex-Biased Genes in Preadult Flies 
Despite the lack of conspicuous phenotypic sexual dimor- 
phism in preadult flies, we found that gene expression differ- 
ences between the sexes were vast compared with the 
differences between developmental stages (fig. 3 and supple- 
mentary table SI and fig. SI, Supplementary Material online). 
The frequency of extremely male-biased genes was higher in 
larvae than in pre-pupae (fig. 1 and supplementary table SI, 
Supplementary Material online), which may reflect a greater 
discrepancy between the sexes in gametogenesis in larvae 
before the onset of ovarian development. In contrast to our 
results, a previous microarray study reported more male- than 
female-biased expression in third instar larvae and did not 
detect decreased X-linked male-biased expression (Meisel 
et al. 2012a), potentially due to the differences between 
whole-body and tissue-specific microarray and RNA-seq anal- 
yses. Lebo et al. (2009) took a different approach by contrast- 
ing expression between wild-type flies and mutants lacking 
germline tissue, using microarrays, and report higher levels of 
sex-biased expression in pupae that are closer to what we 
observed in larvae and pre-pupae. 

Continuity and Specificity in Sex Bias 
In contrast to ontogenetic studies of sex bias in birds (Mank 
et al. 2010), where sex-biased expression is largely discontin- 
uous among developmental stages, the D. melanogaster 
genome shows a complex mixture of genes with both con- 
tinuous and stage-specific sex-biased expression; in fact, 
roughly 50% of genes fell into each group (fig. 3). Genes 
were more likely to retain male-biased than female-biased 
expression throughout development, which probably reflects 
the earlier onset of male gametogenesis (Bodenstein 1950). 
Consistent with this interpretation, continuously male-biased 
genes were enriched for GO terms associated with gameto- 
genesis, whereas larval- and pre-pupal-specific male-biased 
genes were not enriched for any terms (supplemental table 
S6, Supplementary Material online). 

Our finding that more genes switched from female to male 
expression than vice versa is consistent with the pattern in 
silkworm (Zhao et al. 2011). We found that adult male flies 
appeared to co-opt genes encoding extracellular components 
in larval females. Further functional characterization of these 



genes is needed to understand both the mechanistic basis of 
reversals in sex bias and the conservation of this pattern 
across taxa. 

In contrast to our results. Can et al. (2010) found higher 
stage specificity of male-biased genes compared with female- 
biased genes, despite the broad consistency between our data 
and theirs noted above. The difference might arise because 
the wild-type flies used in their study were less than 1 day 
posteclosion, such that oogenesis may not have been fully 
activated. 

Evolutionary Divergence in Ontogenetic Context 
The pattern of evolutionary divergence for sex-biased genes 
changed considerably when we accounted for developmental 
patterns. When we ignored the consistency of expression 
across developmental stages and considered divergence 
rates of larval and pre-pupal sex-biased genes, we obtained 
similar results to those reported for adult sex-biased genes 
(reviewed by Ellegren and Parsch 2007; Parsch and Ellegren 
2013); namely, more rapid divergence for male-biased genes 
through relaxed purifying selection, slower divergence for 
female-biased genes with weak evidence for stronger purifying 
selection, and a fast-X effect (table 1 ). 

However, when we incorporated developmental patterns, 
our analysis of genes with continuous or stage-specific sex 
bias in expression shows that patterns of evolutionary diver- 
gence of sex-biased genes are more nuanced than previously 
appreciated. Our results show that the well-known pattern of 
more rapid evolution of male-biased compared with female- 
biased and unbiased genes does not always hold when the 
developmental pattern is considered (fig. 4). Furthermore, 
male- and female-biased genes can evolve at similarly rapid 
rates, depending on their developmental pattern, but that 
stage limitation influences divergence rates in distinct ways 
in males and females (fig. 4b and d). This finding suggests that 
the sexes exhibit distinct patterns in how selection acts 
throughout development. 

Our finding that stage-specific genes evolved more rapidly 
than continuously expressed genes within each sex supports 
the notion of developmental evolutionary constraints 
imposed by shared gene expression across multiple time 
points (Artieri et al. 2009; Mank et al. 2010), analogous to 
the constraints thought to be imposed by expression across 
multiple tissues (Haerty et al. 2007; Mank et al. 2008; Meisel 
201 1 ). Three further features of these results are notable. First, 
and surprisingly, the difference in divergence rates between 
sex-specific genes disappears when considering the continuity 
of expression across stages. Male- and female-expressed genes 
specific to larvae evolved equally rapidly, and genes continu- 
ously expressed in either sex also evolved at similar rates 
(fig. 4). Second, the more rapid divergence of female larval- 
specific genes was consistent with relaxed purifying selection 
acting on these genes, relative to genes continuously 
expressed in females. However, male larval-specific genes, 
though diverging relatively rapidly, did not show this pattern, 
suggesting that distinct evolutionary processes can shape 
the evolution of sex-specific gene expression even when 
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divergence rates are similar. Finally, genes specific to third 
instar larvae within both sexes evolved most rapidly, with 
genes specific to adults evolving relatively slowly. This result 
contrasts with Artieri et al.'s (2009) finding that adult-specific 
genes evolved more rapidly than a combination of all larvae- 
and pupae-specific genes. Our study samples two precise 
developmental time points using a higher resolution assay 
of gene expression (RNA-seq as opposed to EST libraries), 
and we show that there is in fact not a continuous increase 
in rates of evolutionary divergence throughout development. 

Conclusion 

Overall, our results highlight that new aspects of the evolu- 
tionary dynamics of sex-biased genes and sex-specific selec- 
tion can be uncovered by considering ontogenetic patterns of 
gene expression. Important next steps for progress in these 
areas include measuring the developmental patterns of sex- 
specific isoforms, which can be as frequent as sex-biased genes 
in adults (Can et al. 201 0; Chang et al. 201 1; Daines et al. 201 1 ), 
and characterizing early somatic sex differences in expression, 
which can be considerable in adults (Chang et al. 2011). 

Materials and Methods 

Experimental Animals and Developmental Stages 
A fly stock was obtained in which the female-specific pro- 
moter Sxl was linked to a gene expressing GFP, such that 
female embryos express GFP and embryos can be sorted by 
sex (Graham et al. 2006). This construct was backcrossed into 
a genetically variable fly population (LHm, Chippindale et al. 
2001; backcrossing performed by U. Friberg). Flies were main- 
tained at 25 °C and a 12:12 lightdark schedule. 

Gonad imaginal discs or (in adults) mature gonads were 
dissected from 1) wandering third instar larvae that had dark 
blue food visible in their guts (i.e., having recently stopped 
feeding and >12 h from pupation, Karim and Thummel 1992; 
equivalent to puff stage 1, Andres and Thummel 1994), 2) 
newly formed white prepupae, and 3) 24-h post-eclosion 
virgin adults. 

Sample Preparation 

The parents of flies used in this experiment were grown at a 
standardized low density (25/ml food) to minimize variation 
associated with parental rearing conditions. These parental 
flies were 3-5 days postemergence when we collected eggs for 
the flies to be sampled. The larvae and pre-pupae we sampled 
were grown at low density (25/ml food) on standard food 
media containing bromophenol blue (0.05%) to facilitate 
staging (Maroni and Stamey 1983). Dissections were per- 
formed from 2 to 7h Zeitgeber time. Animals were sexed 
by fluorescence under a fluorescence microscope and by 
visual inspection of pre-gonads (Kerkis 1931). The gonadal 
discs were dissected on ice-cold phosphate buffered saline 
and immediately transferred to QIAzol lysis reagent 
(Qiagen, West Sussex, United Kingdom) on ice, followed by 
tissue rupture and storage at — 80°C. Although gonad tissue 
is comprised of somatic and germline cells, germline cells 
contribute the majority of mRNA (Gupta et al. 2006). 



We sequenced three nonoverlapping pools (i.e., biological 
replicates) of flies per sex for both larvae and pre-pupae, along 
with one pool for each sex for adults to enable us to compare 
our results with existing data sets. We included more female 
than male pre-gonads within each pool for larvae and pre- 
pupae to reflect the higher RNA yield from male pre-gonads. 
The number of pairs of gonads ranged from 39 to 104 for 
larval and prepupal pools and was 10 for adult pools. 

RNA was extracted using a lipid tissue kit (Qiagen RNeasy, 
West Sussex, United Kingdom), following manufacturer's 
instructions and including the DNase digestion step. mRNA 
was prepared for sequencing using standard lllumina 
protocols. 

Sequencing and Bioinformatics 

Library preparation and RNA sequencing (lllumina HiSeq 
2000) was performed by the Wellcome Trust Centre for 
Human Genetics, Oxford. Each biological replicate was bar- 
coded and pooled on a single lane, yielding 37.17 Gb of 100 bp 
paired-end reads. Read quality and potential sequencing 
biases were assessed using FastQC (http://www.bioinformat- 
ics.bbsrc.ac.uk/projects/fastqc, last accessed August 31, 2013). 
Trimmomatic (Lohse et al. 2012) was used to remove any 
residual lllumina adaptor sequence, trim leading, and trailing 
bases with a Phred quality score less than 4, and using a sliding 
window approach trim when the average Phred quality score 
over four bases dropped below 15. Any paired-end reads for 
which either read was shorter than 25 bp after trimming were 
removed from the data set. After filtering, there were on 
average 13 million reads per pool (biological replicate). 

The genome of Drosophila melanogaster (Adams et al. 
2000) was obtained from EnsembI (release 64; Flicek et al. 

2012) . Each pool's filtered paired-end reads were mapped 
to the reference genome using RSEAA (Li and Dewey 2011), 
which leverages the short-read aligner bowtie (Langmead 
et al. 2009). rRNA genes and rRNA pseudogenes were re- 
moved, as rRNA can bias expression estimates of mRNA 
genes. On average, 81.11% of reads mapped to the reference 
genome per pool and a total of 14,705 genes had mapped 
reads. A filter of four reads per million mappable reads 
(pmmr) was applied to remove nonexpressed and lowly ex- 
pressed genes (Harrison, Mank et al. 2012; Moghadam et al. 

2013) . Genes were filtered if they did not pass the four-pmmr 
threshold in at least two out of three pools for either sex of 
either larvae or pupae pools, or did not pass the threshold 
for either adult pool, resulting in 10,587 genes that were ex- 
pressed in at least one sex and life stage. The raw read count 
outputs from RSEAA were quantile normalized. Chromosome 
location and GO terms (Ashburner et al. 2000) were obtained 
for each gene, if available, from the FlyBase gene associations 
(release FB2012_01; McQuilton et al. 2012). 

Gene orthology information was obtained from FlyBase 
and one-to-one orthologs identified for the melanogaster 
subgroup (Drosophila 12 Genomes Consortium 2007; 
Larracuente et al. 2008), which includes D. melanogaster, 
D. sechellia, D. simulans, D. yakuba, D. erecta, and D. ananas- 
sae. Within the melanogaster subgroup, 9,162 orthologs were 
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identified in FlyBase, and each orthologous group of six genes 
was aligned with PRANK (Loytynoja and Goldman 2005). 
Nine of these orthologous groups were removed from the 
data set due to a member's codons not being in multiples of 
three. A further 27 orthologous groups had erroneous FlyBase 
entries and were missing one or more of the melanogaster 
group species genes. Repeatmasker (http://www.repeatmas 
ker.org, last accessed August 31, 2013) identified 882 genes 
as containing retrotransposons, including long interspersed 
nuclear elements (LINEs) and long terminal repeats (LTRs), or 
DNA transposons, and these groups were removed from the 
data set. In-frame stop codons were found in 124 orthologous 
groups genes and a further 13 orthologous groups were ex- 
cluded because genes were less than 100 bp in aligned gapless 
length. These filtering steps left 8,1 12 one-to-one orthologs of 
the six melanogaster group species. 

These 8,1 12 genes were analyzed with PAML (version 4.4b; 
Yang 2007), using the melanogaster subgroup phylogeny 
(Drosophila 12 Genomes Consortium 2007; Larracuente 
et al. 2008). Along the D. melanogaster branch, the number 
of nonsynonymous substitutions (Dp), the number of poten- 
tial nonsynonymous sites (N), the number of nonsynon- 
ymous substitutions per nonsynonymous site (d^), the 
number of synonymous substitutions (Ds), the number of 
potential synonymous sites (S), the number of synonymous 
substitutions per synonymous site (ds), and their ratio {d^/d^ 
were extracted for each orthologous group. These values were 
verified against d^/ds results from several previous studies 
(Parisi et al. 2003; Ranz et al. 2003; Gibson et al. 2004; Parisi 
et al. 2004; Stole et al. 2004; Mclntyre et al. 2006; Goldman and 
Arbeitman 2007; Wayne et al. 2007; Mikhaylova et al. 2008; 
Ayroles et al. 2009; Can et al. 2010; Innocenti and Morrow 
2010; Wasbrough et al. 2010; Wyman et al. 2010; Graveley 
et al. 2011) obtained from the sex bias database SEBIDA (ver- 
sion 3.0; Gnad and Parsch 2006). 

Data on standing variation in protein-coding sequences 
were obtained from the Drosophila Genome Reference 
Panel (DGRP; Database file date July 18, 2012; Mackay et al. 
2012) to calculate the numbers of synonymous (Ps) and non- 
synonymous polymorphisms (Pp). For the 8,112 genes used 
for the divi/ds analysis, SNPs were extracted from the DGRP 
database coded as either synonymous or nonsynonymous, 
with an additional filter requiring the minor allele to be pre- 
sent in at least 10% of lines (17/168). These data were used to 
calculate the direction of selection statistic, DoS = Dn/ 
(Dp + Ds) - Pn/(Pn + Ps) (Stoletzki and Eyre- Walker 2011). 

Analyses 

Within a developmental stage, we defined sex-biased genes as 
having at least 2-fold differential expression and sex-limited 
genes as having zero expression in one sex. For both catego- 
ries, we required a significant P value from a t-test after mul- 
tiple testing correction (P^dj < 0.05, Benjamini Hochberg 
correction [Benjamini and Hochberg 1995], false discovery 
rate [FDR] 5%). This approach to categorizing sex-biased 
genes was consistent with two alternative approaches 
(EdgeR and DESeq): our categorization showed an 87-95% 



overlap with these methods, which compares well with the 
overlap observed between them (79-93%). We also explored 
patterns of sex-biased gene expression using different thresh- 
old of differential expression. 

To examine the continuity of sex-biased gene expression 
throughout development, we used two methods to catego- 
rize adult genes. First, we used our own adult data, following 
the above procedure, requiring P^dj < 0.05 in a genome-wide 

test. Second, we used categories defined by the SEBIDA 
meta-analysis with a 2-fold threshold (Gnad and Parsch 2006). 

In all analyses, we treated autosomes and X chromosomes 
separately because they display different patterns of sex- 
biased expression and rates of evolutionary divergence in 
D. melanogaster (e.g., Meisel et al. 2012a). 

To test for differences in d^/ds among groups of genes, we 
used permutation tests with 1,000 iterations. To test for dif- 
ferences in DoS, we used Wilcoxon nonparametric tests. We 
present bootstrapped confidence intervals (based on 1,000 
replicates) with d^lds and DoS values, t-tests, correlations, 
tests, and Wilcoxon tests were calculated using JMP v. 10 
(SAS Institute Inc., Gary, NC, USA). Permutation tests and 
bootstrapping were conducted in Excel 2010 (Microsoft, 
Redmond, Washington, USA). 

We used GOrilla (Eden et al. 2009) to test for GO term 
enrichment in groups of sex-biased genes. 

Supplementat7 Material 

Supplementary analyses, supplementary figures S1-S4, and 
supplementary tables S1-S8 are available at Molecular 
Biology and Evolution online (http://www.mbe.oxfordjour 
nals.org/). 
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